A DUAL ALGORITHM FOR NON-ABELIAN YANG-MILLS COUPLED TO 

DYNAMICAL FERMIONS 



J. WADE CHERRINGTONl 



00 
O 
O 

(N 



m 
(N 



o 



^Department of Applied Mathematics, University of Western Ontario, London, Ontario, Canada 

e-mail: jcherrin@uwo.ca 



Abstract. We extend the dual algorithm recently described for pure, non-abelian Yang-Mills on the lattice 
to the case of lattice fermions coupled to Yang-Mills, by constructing an ergodic Metropolis algorithm for 
dynamic fermions that is local, exact, and built from gauge-invariant boson-fermion coupled configurations. 
^ ■ For concreteness, we present in detail the case of three dimensions, for the group SU (2) and staggered 

O ' fermions, however the algorithm readily generalizes with regard to group and dimension. The treatment of 

I the fermion determinant makes use of a polymer expansion; as with previous proposals making use of the 

polymer expansion in higher than two dimensions, the critical question for practical applications is whether 
the presence of negative amplitudes can be managed in the continuum limit. 

1. Background 

. 

Despite continued progress in algorithms and hardware, the inclusion of dynamical fermions in lattice 
Q^' gauge calculations continues to incur significant computational expense. To motivate our proposal for a 

O I novel fermion algorithm, we briefly review how dynamical fermions are currently addressed. Recall that 

dynamic fermions coupled to a gauge field on a D-dimensional hypercubic lattice for > 2 are governed by 
an action of the form 

(1) S[ge,^v,i'v] = SG[ge]+ SF[ge,i'v,i'v] 



where the ge are valued in the gauge group G at the edges of the lattice and ipv are the fermion fields defined 
at the vertices of the lattice. 



■ Unlike gauge group variables, it is not practical to directly simulate Grassmann variables on the com- 



puter. A common approach to dynamical fermion simulation starts by integrating out the fermion variables 
appearing in Splgej^vTipv], to give a function of the gauge variables known as the fermion determinant (its 
I specific form is reviewed in Section 2). The fermion determinant can be combined with the kinetic part 

■ of the gauge boson amplitude e~^'^^^^^ to give an effective action for the gauge variables from which sim- 
ulations on a computer can in principle proceed. However, the fermion determinant renders this effective 

, . ■ action non-local — it couples together gauge variables that are arbitrarily distant in the lattice. This poses 

r> I a considerable problem for the simulation, since computing the change in the effective action due to a small 

■ change in any variable becomes very expensive, growing prohibitively with increasing lattice volume. A 
variety of algorithms have been devised to work with the fermion determinant; a description of some of the 
methods commonly employed can be found for example in [6]. 

After reviewing the description of single component, staggered free fermions in terms of self- avoiding 
polymers (as was done for example in [9j), we review what happens when a similar procedure is applied to 
multi-component fermion fields minimally coupled to gauge fields. In this case each polymer configuration 
corresponds to a Wilson loop functional; i.e. the trace of a product of representation matrices around the 
polymer. Because there is more than one component of the fermion fields in the non-abelian coupled case, 
the strict self-avoiding constraint of the single component case is weakened; that is, for an n-component 
fermion, up to n directed polymer lines can enter and leave a given vertex. The picture has long been 
known — it is essentially that of a hopping parameter expansion of the fermion determinant, described for 
example in [l^. Unlike many past applications, in the present case no cut-off in the power of the hopping 
parameter or otherwise is applied. Because we seek an exactly dual model, all polymers are included in the 
configurations considered. 

For each polymer diagram that arises in the free case, upon applying the duality transformation for the 
group- valued field the result is a sum of configurations consisting of all closed, branched, colored surfaces (spin 



foams) with open one-dimensional boundaries defined by the polymer diagram. The totality of spin foams 
associated with all polymer diagrams (including the trivial empty polymer) defines the joint configuration 
space. Crucially, local changes to the dual configurations (either polymer or surface structure) lead to local 
changes in the dual amplitude. 

The two theoretical inputs for this construction, a polymer decomposition of the fermion determinant (i.e. 
hopping parameter expansion as described in [T4j) and a dual non-abelian model (e.g. [T^ and references 
therein), have been present in the literature for some time, and as we shall see the construction of the joint 
dual model at the formal level is a rather straightforward synthesis of these constituent models. However, 
unlike (the simplest implementations of) conventional lattice gauge simulations, finding any practical algo- 
rithm for a dual model has proven somewhat non-trivial in the non-abehan case for dimensions greater than 
two. The algorithm proposed here builds upon the dual non-abelian algorithm of [5] that has recently been 
tested in the pure Yang-Mills sector. In addition to pure spin foam moves, we construct a set of moves that 
act on polymer structure and specify the type of vertex amplitudes that arise due to the charges carried by 
the polymer. Currently, an implementation of this algorithm is being tested and will be reported on in a 
forthcoming work. 

For context, it should be noted that a similar picture was present in the work of Aroca et al. [1] and 
Fort [7], which dealt with the abelian case of U{1) and proposed using a Hamiltonian that leads to a different 
Lagrangian formulation, where the ensemble is built from a restricted subset of the configurations that arise 
in the Kogut-Susskind case. In future work, we beheve the non-abelian generalization of [l] may be a very 
interesting alternative to the Kogut-Susskind formulation used here, particularly if the imbalance between 
negative and positive amplitudes and the reduction of species doubhng described in [l] can be carried over 
to the non-abelian case. 

The outline of the paper is as follows. In Section 2, we review the origin of the fermion determinant 
and discuss its expansion in terms of polymers. In Section 3, we briefiy review the dual computational 
framework for non-abelian, pure Yang-Mills theory on the lattice. In Section 4, we show a natural way to to 
combine these frameworks and formulate ergodic moves for the coupled fermion-boson system. In Section 5, 
we offer some conclusions and describe our program for ongoing numerical work based on this algorithm and 
its extensions. Appendix A describes the vertex amplitudes that arise in the joint case, while Appendix B 
expands on Section 2 to describe the trace structure of polymers with multiply occupied vertices. 



In this section, we start from a conventional lattice discretization of free fermions, following [12] in 
essentials and notation. In the usual manner, the fermion determinant is arrived at by exact integration 
of the Grassmann variables. Following |9], the fermion determinant is then expanded into states, each of 
which is represented by a family of disjoint, closed oriented loops, including trivial and degenerate "loops", 
the monomers and dimers, respectively. A typical polymer configuration (in the D = 2 massive case) is 
illustrated We now review explicitly how the fermion determinant and polymer picture come about. Note 
this section is purely for pedagogical purposes and to fix notation to be used later; those familiar with the 
hopping parameter expansion of the fermion determinant can safely skip it. 

2.1. Kogut-Susskind staggered fermions — free case. To illustrate the concept and introduce termi- 
nology, we treat a single species of fermions with no additional indices. We start with the naieve lattice field 
action for free staggered fermions 



2. Polymer description of fermions on the lattice 



(2) 




where are the (Euclidean) Dirac matrices and V is the set of lattice vertices. 
Using the central difference for the partial derivative, this becomes 



(3) 
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Figure 1 shows the fermion state on part of a larger lattice to which periodic boundary conditions are applied. 
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Cycle Monomer Dimer 



Figure 1. Part of a typical configuration in the polymer expansion of the fermion deter- 
minant in two dimensions (massive case). 




where D is the dimensiorH of the hyper-cubic lattice, m is the mass, a is the lattice spacing, and ^ labels 
one of the D directions of the lattice; /t is the unit lattice vector associated to the /ith direction. Following 
[12] . we change to a basis which diagonaHzes the gamma matrices, rewriting the free action as 

(4) 5= ^ |M(V^,V.)-i^E"-M(V^xV'x+A-f 

with ax^i = (— 1)^^^* i-aij,-! fgj. ^ g 2, D] where the Xi are the components of the lattice site four-vector 
and ^ — ma. The edge dependent sign factor Uxfi arises from the chosen diagonaHzation. 

To express the result of integration over the Grassmann variables it is convenient to introduce the quark 
matrix Q, defined in terms of lattice regularized action as 

(5) SF[ge,fpv,ijy] ^ "^yQidelyxtpx- 

For later reference, we include dependence on the gauge degrees of freedom in our definition of Q. We next 
apply the well known result [12] that the Grassmann integral over the fermion fields at every vertex evaluates 
to 

(6) flu d^P^chpA e-^Ag.,^.,-^.] =/"[[] d^'vcbpA ^-«e^^»'3[9e]„.^. ^ detQlg^], 
•' \vev / '' \vev J 

the determinant of the quark matrix. 

Next we recall the continuum form of the action for massive fermions coupled to a gauge field. In the 
massive case, the quark matrix can be written as 

(7) Qyx — ^I^yx Kyx: 



In three dimensions, the continuum limit of the staggered fermion action contains flavours correponding to two inequivalent 
representations of the Dirac algebra [2]. 



where M is the fermion mass and K^y is the hopping matrix that is non-zero for nearest neighbor pairs 
{x,y). By inspection of (U, we write the quark matrix as 

D 

(8) Qyx = M5yx ~ K'^ax^,{5y,x-[i - ^y-^x)- 

tj.=i 

We now apply a well known identity for the determinant of (any) matrix Q, 

(9) det Q = ^ sgn(7r) Y[ Qxtt(x) , 

7T X 

where tt ranges over the set 11 of all permutations of the indices of Q and sgn(7r) is the sign of the permutation. 
In the case of the quark matrix Q, the matrix indices being permuted correspond to vertices of the lattice. 
Thus, one is led to consider the product Yix Qx-k(x) for every permutation of lattice vertices. 

The next step is to recognize that every permutation tt can be decomposed into a composition of disjoint, 
non-trivial cyclic permutations tt'^, tt = Ili'^'i^- ^^O'" ^ matrix with all entries non-zero, these permutations 
may involve sets of vertices that are arbitrarily separated on the lattice. However, the quark matrices that 
arise in practice have a very specific structure (originating in the lattice discretization from nearest neighbor 
approximations of the derivative operator) ; the only non-zero matrix elements consist of nearest-neighbor 
pairs (and in the massive case, on-diagonal). Thus, the non-zero contributions can be analyzed as follows. 

A given permutation tt affects any vertex trivially (the vertex is sent to itself) or as part of a non-trivial 
cyclic permutation. In the massive case, vertices that are permuted trivially give monomer factors equal 
to the mass M . In the massless case where diagonal entries vanish, any trivial permutation will lead to a 
vanishing contribution; thus every vertex must participate in a cyclic permutation in the massless case. 

For the non-trivial permutations, it is useful to distinguish two cases that a vertex may participate 
in. Permutations that swap a pair of neighboring vertices are referred to as dimers; all other non-trivial 
permutations consist of non-trivial loops of edges on the lattice; we shall refer to these as cycles. Given these 
observations on the structure of Q, we can now write an expression for detQ in more explicit detail as 

(10) detQy, = ^sgn(7r)Af^" Y[ Kxy. 

Where Nm are the number of monomers. The product is over all directed edges {xy) that are part of a dimer 
or cycle of the permutation. 

2.2. Coupled case. To couple fermions to the gauge fields, the ordinary derivative is replaced by the 
covariant one, thereby introducing the gauge variables ge which act on the fermions through the matrices 
of the representation corresponding to the charge of the fermion. For specificity, we will consider fermions 
charged in the representation of G labelled by c. One can show [l^ that the lattice action for staggered 
fermions coupled to the gauge field becomes 

(11) ^ = E U^(^x^-) + KY.ax^ (?.c/i^^.+A - i^x+f.Ux^i'.yj ■ 

X \ / 

In contrast to the (simplified one component) free case, there is in general a vector possessing multiple 
Grassmann variable components at each vertex, and matrices U{ge) that act non-trivially on this vector. In 
terms of our permutation expansion, the quark matrix now takes on component as well as vertex labels as 
follows: 

(12) ^j^[5e,^.,?J = E ^y^W-^- 

x.y^V 

Comparing lfT2|) to (fTTj) we identify the multicomponent quark matrix as 

D 

(13) Qifxige] = MSyx -Kj2ax, {{UlrSy,x-^ - iUxf.rSy^f.,x) • 

The ge dependence is through the representation matrices Ue, where e is labelled in one of the two conven- 
tions introduced above. The determinant formula can again be applied; upon doing so, the expansion into 



Figure 2. Graphical representation of one (solid) of all possible permutations (dashed) 
associated with a cycle of vertices of length L. 




Observe that the permutations tt now act on both vertex and component indices of Q; the action of tt 
on indices and vertices can be separated into the maps tt/ and ny, respectively. As in the free case, the 
locality structure allows one to identify non- vanishing, polymer-like contributions. For tt where vertex index 
mapping is trivial, the only non- vanishing entries are those for which 7r/(a;, i) = i as the massive term is an 
inner product with no cross-terms. Thus, the monomers contribute factors Mn, where n are the number of 
components of the fermion vectors. 

As in the single-component case, for vertices that participate in non-trivial permutations, one still has 
that only permutations which move x — > tt{x) where Tr{x) is a nearest neighbor are non- vanishing. However, 
due to the presence of multiple components to the fermion field, permutations can shift both components at 
a vertex simultaneously. 

Configurations which involve non-trivial shifts of more than a single component (multiply occupied vertex 
contributions) are discussed in the Appendix B. In the remainder of this section, we will restrict ourselves 
to the case where only a single component participates in the shift, to illustrate in a simple setting how the 
trace of a product of Ue matrices comes about. 

By inspection of (fTTI) . we see for a given nearest neighbor vertex shift, all possible permutations of indices 
are allowed, since in the general case U^^^ has all non-vanishing entries. To continue our analysis we factor 
the permutation into a part that acts on vertices Try (these correspond to the dimers and cycles of the free 
case) and a part that acts on component indices ttj. We now write the fermion determinant as 

(15) det Q[ge\ = sgn(7rv') TT Q^^^'^^/^} + (multiply occupied vertex contributions) . 



Focusing our attention on the first term, we note that only one component per vertex is shifted in each of the 
products of the sum (multiple component shifts at a vertex are precisely what is included in the second term, 
discussed in Appendix B). For a given Try, we can represent tti as an ordered sequence of arrows through the 
discrete n-point space living above every vertex acted on by Try; see Figure [21 Note that for a given dimer 
or cycle Try, all possible index sequences correspond to permutations of the same order (thus sgn(Try) can be 
factored out). The final step in the analysis is to recognize that the sum over all paths is simply the trace 
of the matrix product of representation matrices around the cycle. 

We define Di as the restriction of det(5[(7e] to permutations involving singly occupied vertices; that is, 
the first term of lfT5|) . Di can be constructed out of loops with a single associated trace as follows: 



where Ne is the total number of edges where vertices are shifted in the permutation. In the second line, the 
Einstein summation convention for repeated indices is used. Observe that {xy) denotes an oriented edge. 



TT/ X 



(16) 





Figure 3. Visualization of 2-component fermion field on 16"^ lattice. Polymers that appear 
to have open ends close on the opposite side of lattice due to periodic boundary conditions. 



Depending on whether {xy) is along or opposing the canonical orientation, one has either a product of C/(a;j,) 
or U(yx) — Uj^xy) ■ "^^^ visualization of a typical polymer configuration on a = 3 lattice (including those 
with multiply occupied vertices) appears in Figure [3] below. 

So far our discussion has been generic with regard to group, dimension, and fermion charge; the only 
major choice has been staggered fermions rather than an alternative lattice discretization. In the remainder 
of this work, we restrict our attention to D — 3, G = SU{2), and rt-component massive fermion fields charged 
with half- integer spin c and minimally coupled to SU{2). 

3. Spin Foam Description of Pure Gauge Theory on the Lattice 

Having described what we shall refer to as the free (no coupling to gauge fields) fermion partition function 
in terms of closed lattice polymers, in this section we briefiy review the dual formulation of pure (no coupling 
to fermions) Yang-Mills, which leads to closed, colored, branched lattice surfaces. We shall see in the next 
section that these pictures naturally combine to give the full interacting partition function in terms of a 
space of coupled configurations. 

It can be shown (see for example [3,, 13], and [5j for detail on G = SU{2) in three dimensions) that starting 
from the lattice discretized action for pure Yang-Mills 

(17) Zb= f l[dg,e-^pe.S(a.)^ 

'' eeE 

one can transform to a spin foam formulation expressing the partition function in terms of dual variables as 
follows: 

(18) = E (e n i8/(».,j.) n ^^(*e,je)-M f n e-t^-(^'+^H2j, + 1) ] . 

Here V , E, and P denotes the vertices, edges, and plaquettes of the lattice, respectively. The summations 
over i and j range over all possible edge and plaquette labellings, respectively. A plaquette labelling j 
assigns an irreducible representation of SU (2) to each element of P. These representations are labelled by 
non-negative half-integers (we will denote this set by ^N), also referred to as spins; a labelling j is thus a 
map j: P ^ ^N. In the SU{2), D = 3 case, edges are also labeled by half-integer representations, thus an 
edge labelling is a map i: E ^ ^N. 

6 



Figure 4. Visualization of pure Yang-Mills vacuum on 16'^ lattice. 



Following [5], we define a (vacuum) spin foam configuration as one summand in (fT8| . i.e. a labelling of 
both plaquettes and edges by spins and intertwiners, respectively. The amplitude assigned to a spin foam 
factors into a local product of amplitudes. The vertex amplitude (18j symbol) depends on the 12 plaquettes 
and 6 edges incident to a vertex; the factors depend on the 4 plaquettes and the intertwiner labelling of 
an edge, and there is a product of local plaquette factors. 

The locality of the spin foam formulation was applied in [5] to perform computations that were verified 
against conventional methods. 



4. Dual Fermion-Boson Simulations 

In this section we describe how the dual pictures of lattice fermions and gauge bosons presented in the 
previous two sections can be combined to form a joint dual partition function, built up of gauge-invariant 
configurations with discrete occupancy and representation labels. 

4.1. The Joint Partition Function. Using the action S[ge, tpv,tpv] for the full theory we write the partition 
function as follows: 



(19) Zj = J (Hdgej J ( n d^^dijA e-^-[9' 



,V'^>,'^«lg-S'G[Se] 



eeE / \vev 
[]dgJdetQ[ge]e-^«l^= 



KeeE 



where we have integrated out the fermionic variables to get the fermion determinant. We next use the 
polymer expansion for the determinant in the case of gauge-coupled Q, as described in Section 

(20) Zj - J (^]ldg}j det Q[ge]e-'^'^^^^^ 

\eeE J 7 \(xy)e-i J \{xy)e', ) 



Here Nk is the number of K factors (one per unit of edge occupancy) in the polymer configuration; the 
definition of polymer configuration for £) = 3, G = SU{2) is given in Section 4.3.1. In the second line, we 



have substituted the polymer expansion for fermion determinant. Next, we recall the form of the character 
expansion (see [5] and references therein) for the amplitude based on the heat kernel action at a plaquette p, 

(21) = L_^(2j- + l)e-^J-0-+i)^^.(^), j = 0,i,l,... 

K{I,\) j 

Substituting the character expansion into the previous equation, we have 

(22) Zj = f mdgAY.sgnmnMf-K''- I J] Tr [ [] f/(.,) 

VeG-E / 7 \(xy)ej J \(xy)G7 

X nE(2jp + l)e-*^'^'^'''+'^X,(5), 

2 

where an overall constant factor of K{I, ^) per plaquette has been discarded. We show in Appendix A that 
the group integrals over products of traces and characters in each term of the character expansion can be 
evaluated exactly in terms of charged 18j symbols, provided a sum over intertwiner labels is made at each 
edge. Using the vertex and edge amplitudes of Appendix A, we can exhibit the joint dual partition function 

as 

(23) zj = EEE(^(T')nw^(*-^-^)n^(*-j-^)"'n^'^'''^''^'^(2j,+i) 

yev j i \ vev eS-E peP 

where 5(7) = sgn(7) J^j^^j^-jg^ a(j:j,)(nM)^"X^^' combines the two sign factors and a product of M and K 
factors. As we shall see in the next section, joint configurations associated to a polymer 7 carry in general 
three rather than a single intertwiner label ie for each edge belonging to the polymer; i here ranges over all 
the intertwiner labels. 

Although the overall dual amplitude is still a product of local amplitudes as in the pure Yang-Mills case, 
the presence of the Wilson loop functionals associated to non-trivial polymers requires the vacuum vertex 
and edge amplitudes to be modified in a way that we define in Appendix A; the result is a product of 
modified 18j^ symbols and edge amplitudes N"^ that are charged according to the polymer content 7 of the 
configuration. 

The joint ensemble that results here can be viewed as a generalization of the usual definition of spin foams 
to include one-dimensional structure corresponding to the presence of fermionic charge. For a given polymer, 
there is a sum over all spin foams satisfying admissibility, which is modified at the polymer edges. Prom the 
worldsheet point of view [4|, a polymer loop acts as the source or sink of 2c fundamental sheets, where c is 
the half-integer charge of the fermion. 

4.2. The Joint Fermion-Boson Configurations. In this section we define explicitly the set of configu- 
rations that include all those that give non-zero contribution^ to the joint dual partition function. 

Specifically, for a given polymer 7, we introduce definitions that will allow us to characterize the set of 
spin foam configurations (plaquette colorings) that are admissible in the presence of 7. 

As in the pure Yang-Mills case [5], we assume a splitting has been made for each edge with ji and j2 on 
one side and and ji on the other. Because of the presence of charges ci and C2, the intertwiner is generally 
6-valent and three splittings have to be made. For discussing edge admissibility, we assume the splitting is 
such that ci and C2 are in the middle of the channel as shown in Figure [H Less symmetric splittings are 
possible, but we restrict our attention to this case in the following. 

Let 7 denote a polymer configuration of charge c. We define the set of admissible plaquette configurations 
as those configurations whose labellings satisfy c-edge admissibility at every edge in the lattice, where c- 
admissibility is defined as follows: 

Definition 4.1 (c-edge admissibility). The spins assigned to plaquettes incident to an edge are said to be 
c-edge admissible if the parity and triangle inequality conditions are satisfied. The charge insertions ci and 



''Due to exceptional zeros there may be configurations that are admissible by the conditions defined in this section, but are 
nonetheless zero. As in the pure Yang-Mills case [5], we assume exceptional zeros are sufficiently isolated that ergodicity on 
admissibles is equivalent to ergodicity on non-zero configurations. 



'2 



J 3 U 

Figure 5. Symmetric splitting of a 6-valent SU{2) spin network vertex. 

C2 may be or c, according to whether the edge is uncharged, singly charged, or doubly charged. Writing 
ji; J2, ja and for the four spins incident to a given edge, these conditions are 

(1) Parity: 

3i + 32 + 33 + 34 + C1 + C2 is an integer. 

(2) Triangle Inequality: for each permutation x = {xi,X2,X3,Xi,X5,xe) of the charge and spin vari- 
ables (ci, C2, ji, j2, js, j4) we have 

Xi + X2 + X3 + X4 + X5 > Xq. 

These conditions are equivalent to the existence of a non-zero invariant vector in the SU{2) representation 

ci (8) C2 (8) ji ® J2 <8) js (8> ji- 

The allowed range of intertwiner labels ii, ic and i2 depend on the incident spin labels, on each other and 
on c through the vertices where the charges ci and C2 enter the diagram. We now state the condition for and 
admissible spin foam (plaqeutte and edge intertwiner labelling) in the presence of an arbitrary polymer 7: 

Definition 4.2 (7-admissible spin foam). A spin foam is ^ -admissible if and only if for every edge e € E: 

(1) The plaquettes incident to e are c-edge admissible in the sense of Definition 4.1, with ci and C2 
assigned depending on the occupancy of e by 7. 

(2) Each vertex of Figure[5]is admissible. Explicitly, the following conditions are simultaneously satisfied: 

*2 + ji + j2, ii + ic + ci,i2 + ic + C2 and ii + jz + ji are integers 

and 

ii G - 32\,3i +32] n [\ic - ci|, |ic + ci|], 
ic e [\ii - ci|,ii + ci] n [1^2 - C2I, 1^2 + C2I], 

«2 e [|j3 - 3i\,h +34] n [\ic - C2I, Kc + C2I]. 

For a given polymer 7, we denote the set of 7-admissible spin foams by T^. 

4.3. The Joint Moves and Algorithm. In this section we define moves that transform from one joint 
configuration to another. Together, they are ergodic and obey detailed balance and can thus be used in a 
Metropolis or other Markov chain Monte Carlo algorithms. 

Definition 4.3 (Pure spin foam move). A pure spin foam move consists of a single application of the cube, 
edge, or homology move. In terms of their effect on plaquette spins, these moves are as defined in [5]. 
However, their effect on intertwiner labels needs to be generalized to account for the extra intertwiner 
labellings introduced by polymers. 

Because polymer moves require simultaneous changes in spin foam structure, we use the term "pure" to 
distinguish spin moves that leave the polymer structure unchanged. We now describe the generalization of 
each pure spin foam move to account for extra intertwiner labels. 

Definition 4.4 (Generalized cube move). As described in Definition 2.4 and Appendix A. 3 of [5j. With 
reference to compatible intertwiner moves of Type A, no changes in intertwiner labels are necessary. For 
Type B edges, all three labels are increased or decreased by the same half-unit of spin. 
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Definition 4.5 (Generalized edge move). As described in Definition 2.5 of [5]. Rather than changing the 
single intertwiner label, the three intertwiner labels «i(e), 12(e) and ide) are each randomly changed by 
—2,0, or 2 (to preserve parity) units of spin. 

Definition 4.6 (Generalized homology move). As described in Definition 2.6 of [5], but all three intertwiner 
labels are increased or decreased by one half-unit of spin. 

4.3.1. Polymer moves. In Section [221 we saw how sums over permutations in 11 can be encoded into traces 
of matrix products, ordered according to the orientation of the permutation. Our example was restricted to 
singly occupied vertices, and is generalized in Appendix B. Combining these cases, we see that the sum over all 
permutations in 11 can be represented by traces of products of matrices, if we include diagrams corresponding 
to all possible routings at multiply occupied vertices. This new set of objects, oriented diagrams with routings 
at multiply occupied vertices, we refer to as polymers, and denote by V. It is important to distinguish the 
polymers from their finer-grained constituents, the permutations 11, as polymers can be coupled naturally 
into the spin foam partition function, whereas individual permutations cannot be using the methods here. 

In this section we present a set of moves that are ergodic on the space of polymers T' on a 3-dimensional 
hypercubic lattice. The polymer states at an edge in the 2-component case considered here are as follows. 
Assuming a global orientation has been selected for the edges, an edge can be unoccupied, singly occupied, 
or doubly occupied with the occupied cases carrying both positive and negative orientation. The occupancy 
data at each edge can thus be assigned from the set {—2, —1, 0, 1, 2}. We shall use the term "line of fiux" 
interchangeably with directed polymer line. 

Definition 4.7 (Plaquette move). A plaquette and plaquette orientation (clockwise or counterclockwise) is 
randomly selected. To each edge, a delta occupancy of +1 or —1 (with signs given according to plaquette 
orientation) is assigned, and added to the present occupancy. If the resulting occupancy on any edge 
has magnitude greater than 2, the move is immediately rejected. At each multi-valent vertex, a choice of 
routing is made with equal probability. A proposed move that removes occupancy of edges incident on 
multiply occupied vertices must make further random choice of routing that matches the routing present 
or be rejected. This is necessary to preserve detailed balance. If the move is not rejected, the spin of the 
selected plaquette is randomly decreased or increased by c to satisfy parity. Because a change in the polymer 
occupancy forces a change in both the plaquette spin and the charge structure at an edge, the affected 
intertwiner labels (at each edge of the affected plaquette) must change in a way that is compatible; if not 
the result will be immediately c-edge inadmissible for the new charging. 

This is the most fundamental polymer-changing move, and connects a very large region of the space 
of polymers contributing to the fermion determinant. One can see trivially that the plaquette move can 
create fundamental loops of either orientation when applied to "empty" space (plaquettes of zero occupancy 
edges). Figure [6] illustrates how the plaquette move can deform an existing cycle. The same plaquette move 
of opposite orientation would lead to one of multiple routings with a doubly occupied edge, as shown by 
Figure [To] in Appendix B. 

Because an equally weighted routing choice is made amongst several alternatives when the plaquette move 
of (for example) Figure [TOl is made, the move that reverses a particular routing to get the initial loop back 
should occur with proportional probability, in order to satisfy detailed balance. This will unfortunately lead 
to a lowered acceptance rate, particularly in regimes with high occupancy. 

To identify intertwiner moves compatible with a plaquette move, one must give the choice of splitting 
(grouping of ji, . . . ,j4 into pairs) around the edges of a plaquette, in the same way that intertwiner moves 
compatible with a cube move depends on the splitting (see A. 3 of [5]). For conciseness and generality, we 
have made the present definition of the algorithm splitting independent. Forthcoming work on numerical 
simulations with this algorithm will evaluate alternative splittings and describe the appropriate compatible 
moves. 

Definition 4.8 (Global circle move). For integer valued charge c, a global circle moves adds a single line of 
charge and a minimal cylinder of | charged plaquette spins spanning the lattice. The origin and orientation 
of the cylinder is randomly selected to lie on a plaquette of one of three orthogonal lattice planes through 
the origin. 
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Figure 6. Stretching of a loop by a plaquette move 



For half-integer c, two lines of charge spanning the lattice are added and a minimal surface consisting of 
a line of c-charged plaquettes is added between the charge lines. The position and orientation of the sheet 
is randomly selected to lie on an edge of one of the three orthogonal lattice planes through the origin. 

As in the pure gauge theory case, the non-trivial global topology of lattices with periodic boundary 
conditions leads one to moves that create and destroy structure on a global scale, in this case lines of charge. 

In the half-unit charge case, a second line needs to be added to absorb the flux introduced by the first; 
the "smallest" possible global move places a second Hne of half-unit charge immediately beside the first. 

In the case of integer valued charge, the lowest energy (and hence smallest change in amplitude) structure 
satisfying admissibility is formed by wrapping the smallest possible cylinder supported by the lattice; because 
there are 2n fundamental irreps, n can be incident on one side normal to the line, wrap into a cylinder, and 
enter the other to satisfy c-edge admissibility. 

Definition 4.9 (Dimer move). An edge is randomly selected and a single unit of occupancy is added or 
removed. If the result is not consistent with the presence of dimers and cycle edges (i.e. a cycle that ran 
through the edge is broken by the removal of occupancy), the move is rejected. 

Singly and doubly charged dimers cannot be constructed by plaquette moves. Note a singly charged dimer 
contributes Tv{UeUl) = n with weight while a doubly charged dimer contributes Tv{UeUl)Tv{UeUl) — 
TriJJeUlUeUl) = nP — n with weight (unHke general polymers, we combine the two routings into one 
configuration). Both types of dimers evaluate to constants with respect to the gauge variables, and thus 
don't couple to the gauge bosons. 

Definition 4.10 (Junction move). A vertex is randomly selected, and if multiply occupied, the routing of 
flux is changed. 

A multiply occupied vertex in the case of a 2-component fermion field has two flux paths, which can be 
routed in two different ways. When a multiply occupied vertex first appears as a result of a polymer move, 
one routing is randomly selected (similarly in the inverse case, with weightings to preserve detailed V)alancc 
as discussed above). Thus, the junction moves are stricly speaking unnecessary for ergodicity, but may be 
used to improve performance of Metropolis algorithm. 

4.3.2. The Algorithm. Combining the polymer and spin foam moves, given, we give a statement for a Me- 
tropoHs algorithm ergodic on the joint ensemble. 

Algorithm 4.11. (Joint Fermion-Boson Algorithm). An iteration of the joint algorithm consists of choosing 
one of the seven previously defined moves, which can be organized as as follows: 
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Homology move 



The algorithm can be tuned to improve acceptance rate by adjusting the relative frequency of the attempted 
move types. 

The algorithm also tracks the sign of the configuration changes with the creation and destruction of 
fermion loops, which can occur with any of the polymer moves. The product of monomer factors M and 
hopping factors K are also updated for polymer moves. It is important to emphasize that changes in both 
the sign and other factors require only local consideration of the polymer moves. 

With regard to locality, one sees that the (pure spin foam) homology move and global circle moves will 
lead to updates costing on the order of and L respectively, where L is a characteristic side length of 
the lattice. In the pure Yang-Mills case analyzed numerically in [5], the homology moves have negligible 
influence beyond very small lattice sizes. It remains to be seen how this is modifled in the joint case, and 
how large an influence the global circle moves have. 

Within the scope of the current work, the expectation values of observables depending on dual degrees of 
freedom are computed in the usual manner, by averaging the observable over the Markov chain generated 
by the Metropolis algorithm. For Wilson loop type observables commonly studied, the expectation value is 
actually a ratio of dual charged and dual vacuum partition functions, with static charge corresponding to 
the Wilson loop observable present in the charged partition function and a vacuum partition function given 
by Zj of equation l(23|l . The computation of Wilson loop observables of pure Yang-mills and dynamical 
fermions will be reported on in forthcoming work by the author. 



5. Outlook and Conclusions 

We present here a local, exact algorithm for Metropolis simulation of the fermion-boson vacuum. The 
details have been provided for the case of D = 3 staggered Kogut-Susskind fermions coupled to a Yang- 
Mills SU{2) field; however the algorithm has a straightforward generalization to other dimensions and gauge 
groups. 

A limitation of the algorithm as currently given is the species doubling inherent in the (unrooted) Kogut- 
Susskind formulation (e.g. in four dimensions there will be four species). An alternative to Kogut-Susskind 
fermions which addresses species doubling was developed by Aroca et al. [ij and Fort [7j. We are currently 
investigating this modified fermion action in the non-abelian, higher dimensional context. Another approach 
would be to go through a similar procedure using Wilson fermions, in which the unwanted doublers become 
very heavy in the continuum limit. 

The crucial question for any new fermion algorithm is its performance relative to the highly evolved 
dynamic fermion methods that exist within the conventional lattice gauge community today. In three- 
dimensions, slow-down at weaker coupling has been observed in recent work on the pure Yang-Mills case \5^. 
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The situation in D = 4 is not well understood and is currently the subject of numerical work by the author, 
as are improvements in the original D = 3 case. 

With regard to the continuum limit, we expect the most critical question for the algorithm proposed here 
is the seriousness of the sign problem. A hard sign problem has been discussed as a general feature of the 
polymer expansion in the continuum (small mass) limit [H]. While oscillating signs can be overcome for 
lattice fermions in certain two dimensional theories [TT| [T5] , the author is not aware of methods that have 
successfully addressed the sign problem for D > 2. As both the fermion and dual Yang-Mills (spin foam) 
amplitudes can carry negative signs, an important question is how the signs interact; i.e. the problem of 
signs may be harder or easier than for either the free fermion polymer expansion or dual Yang-Mills alone, 
depending on how the signs correlate. 

Although numerical developments are required to begin evaluating this proposal, we believe the approach 
may be of considerable interest. We find it remarkable that the fermion expansion into polymers and the 
gauge field dualization into spin foams (both of which have been extensively explored on there own) , combine 
together in a way that is very compelling geometrically, and allows a local, exact Metropolis simulation using 
gauge- invariant configurations carrying entirely discrete labels. 

Acknowledgement. The author would like to thank Dan Christensen, Florian Conrady, and Igor Khavkine 
for valuable discussions. The author was supported by NSERC. 



Appendix A. Charged nJ Symbols 



A.l. The dual model with charges. In this appendix, we deal specifically with D — 3, G — SU{2). 
Following the discussion in the appendix of [5], we recall that the dual partition function (in the absence of 
charge) has the form 



(24) 



^^Yi / n ^■^'^ n "^ipxjpidp), 

{jp} ee-B peP 



where summation over jp is over unitary irreducible representations of SU{2). At this point it is convenient 
to specialize to a D = 3 cubic lattice with periodic boundary conditions; orientation choice is as given in 
the appendix of [5]. With this choice of orientation, the holonomy around a plaquette p is gp — gig2g-i^g4^ , 
where 51,52,53 and 54 are the group elements associated to the edges of the plaquette p, starting with an 
appropriate edge and going cyclically. Recall that the inverse g~^ is used if the orientation of edge i does 
not agree with that of p. Thus 

{git U,, (52)g U,^ ig^')i U,^ {g^'Td , 



(25) 



Xjp i9p 



where Uj{g)^ denotes a matrix element with respect to a basis of the j representation. If we insert (|25|) 
into (|24|) and collect together factors depending on the group element 5e, we get a product of independent 
integrals over the group, each of the form 



(26) 



dge U,Aget\ UMil U,,{g-'tl U,,{g-')l 




Here and below we use a graphical notation for tensor contractions, as in [5]. 

Equation (|26| defines a projection operator on the space of linear maps ji (g) js ji ® j2, so it can be 
resolved into a sum over a basis of intertwiners li : ji ® js ji ® j2 



(27) 




\ ^ hi* 
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where the intertwiners I* : ji Cg) j2 ji <Xi are chosen such that the trace (/*/ , 7^) of the composite !*,!{ is 
zero whenever i' ^ i and non-zero if i' — i. The projection property is readily verified. 
We next define Z^, the partition function charged according to the polymer 7, as follows 



(28) 



{jp}-^ eGB \ee7 / peP 



Collecting matrix factors by dependence on edge variable ge , we find in addition to the matrices from the 
four incident plaquettes, a matrix from the edge e with charge c belonging to the polymer 70 



(29) 




As in the pure case, the group integral can be resolved into invariant intertwiners 



(30) 






If, for each edge of the lattice, we fix a term i in the above summation, the intertwiners Jj and /* can 
be contracted with those coming from the other edges, leading to a sum over intertwiner labellings at every 
edge. Observe that at edges occupied by polymers, there is more than a single intertwiner spin label due to 
the additional splittings (see Figure[5]) introduced by the charge lines. At each vertex of the lattice, there will 
be six intertwiners li (some carrying multiple labels), and their contraction can be graphically represented 
as an octahedral network plus additional lines depending on how the polymer passes through the vertex. As 
well, at each edge there will be a normalization factor corresponding to the denominator of equation l|30p . 

We consider first the vacuum case. In this case, each edge carries only a single intertwiner label. The 
result is the ISj symbol central to pure Yang-Mills spin foams. 



(31) 




The vertices are labelled by the directions of the associated lattice edges emanating from the given lattice 
vertex, namely ±a;, ±?/, and ±z. The value of the 18j symbol depends on the choice of basis elements li and 
/*/ in l(30|) . the six summation indices i labelling the edges, and the 12 incident plaquette labels j. 

We now turn to the case where there is a (single) polymer along one or more of the edges incident to a 
vertex. Each charged normalization factor N = {I*,Ii) in the denominator of l(27j) depends on the charge 
c of the fermion at that edge, the intertwiner labels i on that edge, and the labels of the four plaquettes 
incident on that edge. In the presence of external charges, the vacuum 18j is modified depending on whether 



In the case where an edge is doubly occupied, the integral involves a sixth matrix (and the resolving intertwiners an 
additional input and output arrow). It is straightforward to generalize the present analysis to this case; the resulting doubly 
charged I8j symbols are shown in Appendix B. 
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the line of charge proceeds directly through the vertex or turns, leaving in a direction perpendicular to entry 
direction. We call these cases charged 18j symbols and denote them by an overline, 18j(ji,, ie, 7) where the 
additional dependence on gamma reflects how the polymer charge is routed through the octahedral network. 
Typical charged 18 j symbols are shown in Figure El Cases where the direction of the arrow on the charged 




-z -z 



Figure 7. Charged ISj symbols with flux Hnes passing through at right angles (left) and 
straight through (right). 

lines is flipped will also occur, depending on the polymer orientation. Additionally, as discussed in the next 
section, more than a single line of flux can pass through a vertex, leading to charged 18j symbols of the form 
shown in Figures [8l and fTTl 

In implementing numerical code for this algorithm, a choice of splitting (grouping of the four plaquettes 
into (ji, j2) and (^3,^4) pairs on opposite sides of the splitting) is made and each vertex resolved into a 3-valent 
sub-network with up to three non-trivial intertwiner labels, as shown in Figured At this point, recoupHng 
moves (see A. 2 of [5], and references therein) can be used to reduce the spin network to sums and products 
of know spin networks such as the 6j and theta networks, for which eflicient algorithms are available. It 
should be noted, however, that different splittings lead to differing efficiency in implementation, so some care 
and experimentation should be applied to finding an efficient splitting. Specific splitting schemes and their 
performance evaluation will be reported on in forthcoming numerical work by the author and collaborators. 




-z -z 



Figure 8. Charged 18 j symbols with two pairs of flux. Cases with flux not in the same 
plane are also possible. 

Appendix B. Polymers with multiply occupied vertices 

In order to couple to the spin foam representation of the gauge theory, we seek to collect the permutation 
contributions to the fermion determinant into traces of products of Ue matrices around closed, oriented loops 
of edges. The case of permutations where a single component is shifted was discussed above in Section [221 
For polymers where more than one component is shifted at a vertex, recovering a trace formula is somewhat 
more subtle. 

In Figure [9l we illustrate a case where a vertex is multiply occupied; the orientation is such that there 
are two possible routings that resolve the ambiguity at that vertex. Neither diagram by itself corresponds 
to the desired sum of permutation contributions. In the matrix multipHcations and traces (viewed as a sum 
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Figure 9. Two routings associated with a polymer that self intersects once at a point. 



Figure 10. A move introducing a doubly occupied edge, for which there are two distinct routings. 

over all paths around a loop), there are terms in each corresponding to paths that are not permutations. 
However, the same undesired terms occur with opposite sign in the two diagrams (as one involves paths that 
form a single loop, the other paths that lie in two disjoint loops) so the sum of both captures the sum of 
permutations associated with the polymer. 

A similar cancellation occurs when two loops share a single edge, i.e. the edge is multiply occupied. 
Because there are two multiply ocuppied vertices, there are 2^ — A routings possible, however only two are 
topologically distinct; two representatives appear in FigurefTOl The cancellation of unphysical paths between 
the traces over differently routed polymers is well known from the hopping parameter expansion (HPE) of 
the fermion determinant as discussed for example in |i4j. As shown in Figure [TU the presence of a doubly 
charged edge leads to a charged 18j spin network with a 6-valent node. As well, the charged normalization 
factor iV on a doubly charged edge is as given in the denominator of equation l(27l) , but with an additional 
c-charged line parallel to the original c-charged line. 
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Figure 11. A charged 18j symbol containing a 6-valent node incoming from a doubly 
charged edge. 
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